home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / fft / factorize.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-05-05  |  3.6 KB  |  177 lines
  1. /* fft/factorize.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <gsl/gsl_errno.h>
  22. #include <gsl/gsl_fft_complex.h>
  23.  
  24. #include "factorize.h"
  25.  
  26. static int
  27. fft_complex_factorize (const size_t n,
  28.                size_t *nf,
  29.                size_t factors[])
  30. {
  31.   const size_t complex_subtransforms[] =
  32.   {7, 6, 5, 4, 3, 2, 0};
  33.  
  34.   /* other factors can be added here if their transform modules are
  35.      implemented. The end of the list is marked by 0. */
  36.  
  37.   int status = fft_factorize (n, complex_subtransforms, nf, factors);
  38.   return status;
  39. }
  40.  
  41. static int
  42. fft_halfcomplex_factorize (const size_t n,
  43.                    size_t *nf,
  44.                    size_t factors[])
  45. {
  46.   const size_t halfcomplex_subtransforms[] =
  47.   {5, 4, 3, 2, 0};
  48.  
  49.   int status = fft_factorize (n, halfcomplex_subtransforms, nf, factors);
  50.   return status;
  51. }
  52.  
  53. static int
  54. fft_real_factorize (const size_t n,
  55.             size_t *nf,
  56.             size_t factors[])
  57. {
  58.   const size_t real_subtransforms[] =
  59.   {5, 4, 3, 2, 0};
  60.  
  61.   int status = fft_factorize (n, real_subtransforms, nf, factors);
  62.   return status;
  63. }
  64.  
  65.  
  66. static int
  67. fft_factorize (const size_t n,
  68.            const size_t implemented_subtransforms[],
  69.            size_t *n_factors,
  70.            size_t factors[])
  71.  
  72. {
  73.   size_t nf = 0;
  74.   size_t ntest = n;
  75.   size_t factor;
  76.   size_t i = 0;
  77.  
  78.   if (n == 0)
  79.     {
  80.       GSL_ERROR ("length n must be positive integer", GSL_EDOM);
  81.     }
  82.  
  83.   if (n == 1)
  84.     {
  85.       factors[0] = 1;
  86.       *n_factors = 1;
  87.       return 0;
  88.     }
  89.  
  90.   /* deal with the implemented factors first */
  91.  
  92.   while (implemented_subtransforms[i] && ntest != 1)
  93.     {
  94.       factor = implemented_subtransforms[i];
  95.       while ((ntest % factor) == 0)
  96.     {
  97.       ntest = ntest / factor;
  98.       factors[nf] = factor;
  99.       nf++;
  100.     }
  101.       i++;
  102.     }
  103.  
  104.   /* deal with any other even prime factors (there is only one) */
  105.  
  106.   factor = 2;
  107.  
  108.   while ((ntest % factor) == 0 && (ntest != 1))
  109.     {
  110.       ntest = ntest / factor;
  111.       factors[nf] = factor;
  112.       nf++;
  113.     }
  114.  
  115.   /* deal with any other odd prime factors */
  116.  
  117.   factor = 3;
  118.  
  119.   while (ntest != 1)
  120.     {
  121.       while ((ntest % factor) != 0)
  122.     {
  123.       factor += 2;
  124.     }
  125.       ntest = ntest / factor;
  126.       factors[nf] = factor;
  127.       nf++;
  128.     }
  129.  
  130.   /* check that the factorization is correct */
  131.   {
  132.     size_t product = 1;
  133.  
  134.     for (i = 0; i < nf; i++)
  135.       {
  136.     product *= factors[i];
  137.       }
  138.  
  139.     if (product != n)
  140.       {
  141.     GSL_ERROR ("factorization failed", GSL_ESANITY);
  142.       }
  143.   }
  144.  
  145.   *n_factors = nf;
  146.  
  147.   return 0;
  148. }
  149.  
  150.  
  151. static int 
  152. fft_binary_logn (const size_t n)
  153. {
  154.   size_t ntest ;
  155.   size_t binary_logn = 0 ;
  156.   size_t k = 1;
  157.  
  158.   while (k < n)
  159.     {
  160.       k *= 2;
  161.       binary_logn++;
  162.     }
  163.  
  164.   ntest = (1 << binary_logn) ;
  165.  
  166.   if (n != ntest )       
  167.     {
  168.       return -1 ; /* n is not a power of 2 */
  169.     } 
  170.  
  171.   return binary_logn;
  172. }
  173.  
  174.  
  175.  
  176.  
  177.